
module mod_grid  
  implicit none
  real(8),dimension(:,:,:),allocatable:: &
       x,y,z, &
       g11,g12,g13,g22,g23,g33, &
       jac
contains

subroutine mod_grid_init()
  use mod_para
  implicit none
  integer,parameter::cfgunit=100,iunit=101
  integer::nbl,dims(3),ierr,i,j,k
  real(8)::a1x,a1y,a1z,a2x,a2y,a2z,a3x,a3y,a3z,sumcg
  real(8)::ca1x,ca1y,ca1z,ca2x,ca2y,ca2z,ca3x,ca3y,ca3z

  character(len=100) gridfilename

  open(cfgunit,file='fdtd.cfg',action='read')
  read(cfgunit,*)gridfilename
  close(cfgunit)

  open(iunit,file=gridfilename,action='read')
  read(iunit,*)nbl
  read(iunit,*) dims
  nx1=2
  nx2=dims(1)-1
  ny1=2
  ny2=dims(2)-1
  nz1=2
  nz2=dims(3)
  allocate(x(nx1-1:nx2+1,ny1-1:ny2+1,nz1-1:nz2+1),stat=ierr);x=0.0;
  allocate(y(nx1-1:nx2+1,ny1-1:ny2+1,nz1-1:nz2+1),stat=ierr);y=0.0;
  allocate(z(nx1-1:nx2+1,ny1-1:ny2+1,nz1-1:nz2+1),stat=ierr);z=0.0;
  allocate(g11(nx1-1:nx2+1,ny1-1:ny2+1,nz1-1:nz2+1),stat=ierr);g11=0.0;
  allocate(g12(nx1-1:nx2+1,ny1-1:ny2+1,nz1-1:nz2+1),stat=ierr);g12=0.0;
  allocate(g13(nx1-1:nx2+1,ny1-1:ny2+1,nz1-1:nz2+1),stat=ierr);g13=0.0;
  allocate(g22(nx1-1:nx2+1,ny1-1:ny2+1,nz1-1:nz2+1),stat=ierr);g22=0.0;
  allocate(g23(nx1-1:nx2+1,ny1-1:ny2+1,nz1-1:nz2+1),stat=ierr);g23=0.0;
  allocate(g33(nx1-1:nx2+1,ny1-1:ny2+1,nz1-1:nz2+1),stat=ierr);g33=0.0;
  allocate(jac(nx1-1:nx2+1,ny1-1:ny2+1,nz1-1:nz2+1),stat=ierr);jac=0.0;

  read(iunit,*) (((x(i,j,k), i=1,dims(1)), j=1,dims(2)),k=nz1,nz2+1), &
       (((y(i,j,k), i=1,dims(1)), j=1,dims(2)),k=nz1,nz2+1), &
       (((z(i,j,k), i=1,dims(1)), j=1,dims(2)),k=nz1,nz2+1)
  close(iunit)

  if(x(nx1,ny1,nz1)>z(nx1+1,ny1,nz1)) then
  x=-x
  endif

  if(y(nx1,ny1,nz1)>z(nx1,ny1+1,nz1)) then
  y=-y
  endif

  if(z(nx1,ny1,nz1)>z(nx1,ny1,nz1+1)) then
  z=-z
  endif

  x(:,:,nz1-1)=x(:,:,nz1)
  y(:,:,nz1-1)=y(:,:,nz1)
  deltah=minval(z(:,:,nz1+1)-z(:,:,nz1))
  z(:,:,nz1-1)=z(:,:,nz1)-deltah

  delta=huge(1.0)
  do k=nz1,nz2
     do j=ny1,ny2
        do i=nx1,nx2
           a1x=x(i+1,j,k)-x(i,j,k)
           a1y=y(i+1,j,k)-y(i,j,k)
           a1z=z(i+1,j,k)-z(i,j,k)
           a2x=x(i,j+1,k)-x(i,j,k)
           a2y=y(i,j+1,k)-y(i,j,k)
           a2z=z(i,j+1,k)-z(i,j,k)
           a3x=x(i,j,k+1)-x(i,j,k)
           a3y=y(i,j,k+1)-y(i,j,k)
           a3z=z(i,j,k+1)-z(i,j,k)
           g11(i,j,k)=a1x*a1x+a1y*a1y+a1z*a1z
           g12(i,j,k)=a1x*a2x+a1y*a2y+a1z*a2z
           g13(i,j,k)=a1x*a3x+a1y*a3y+a1z*a3z
           g22(i,j,k)=a2x*a2x+a2y*a2y+a2z*a2z
           g23(i,j,k)=a2x*a3x+a2y*a3y+a2z*a3z
           g33(i,j,k)=a3x*a3x+a3y*a3y+a3z*a3z
           jac(i,j,k)=a1z*(-a2y*a3x+a2x*a3y)+a1y*(a2z*a3x-a2x*a3z)+a1x*(-a2z*a3y+a2y*a3z)


           ca1x=a2y*a3z-a2z*a3y 
           ca1y=a2z*a3x-a2x*a3z
           ca1z=a2x*a3y-a2y*a3x
           ca2x=a3y*a1z-a3z*a1y 
           ca2y=a3z*a1x-a3x*a1z
           ca2z=a3x*a1y-a3y*a1x
           ca3x=a1y*a2z-a1z*a2y 
           ca3y=a1z*a2x-a1x*a2z
           ca3z=a1x*a2y-a1y*a2x

           sumcg=(ca1x*ca1x+ca1y*ca1y+ca1z*ca1z &
                +ca2x*ca2x+ca2y*ca2y+ca2z*ca2z &
                +ca3x*ca3x+ca3y*ca3y+ca3z*ca3z &
                +2.0*(ca1x*ca2x+ca1y*ca2y+ca1z*ca2z) &
                +2.0*(ca1x*ca3x+ca1y*ca3y+ca1z*ca3z) &               
                +2.0*(ca2x*ca3x+ca2y*ca3y+ca2z*ca3z) )

           delta=min(delta,jac(i,j,k)/sqrt(sumcg))
        enddo
     enddo
  enddo
  delta=delta*sqrt(3D0)


end subroutine mod_grid_init

end module mod_grid
